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We try to improve the Thomas-Fermi model for the total energy and electron density of atoms and 
molecules by directly modifying the Euler equation for the electron density, which we argue is less 
affected by nonlocal corrections. Here we consider the simplest such modification by adding a linear 
gradient term to the Euler equation. For atoms, the coefficient of the gradient term can be chosen 
so that the correct exponential decay constant far from the nucleus is obtained. This model then 
gives a much improved description of the electron density at smaller distances, yielding in particular 
a finite density at the nucleus that is in good qualitative agreement with exact results. The cusp 
condition differs from the exact value by a factor of two. Values for the total energy of atomic 
systems, obtained by coupling parameter integration of the densities given by the Euler equation, 
are about as accurate as those given by very best Thomas-Fermi- Weizsacker models, and the density 
is much more accurate. Possible connections to orbital-free methods for the kinetic-energy functional 
in density functional theory are discussed. 
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I. INTRODUCTION 

Density functional theory (DFT) has proved to be 
a powerful way to determine equilibrium and dynamic 
properties of atoms, molecules, and extended systems 
Q IE ^ 111 ' The modern conceptual framework was es- 
tablished by Hohenberg and Kohn (HK) , who proved 
that the exact ground state energy of a system of N elec- 
trons is a functional E[p] of the ground state electron 
density p(r) only, where 



E[p]^ / p{r)V,,tiv)dv + F[p], 



(1) 



and F[p] is an unknown but universal functional of the 
density and is independent of the external potential 
K.t(r). 

Kohn and Sham (KS) 0, H showed that F[p\ can be 
usefully partitioned into two dominate (kinetic and po- 
tential energy) terms and a residual correction: 



|r-r'| 



dr'dr + E,,[p]. (2) 



The first term on the right of eq (O is the KS kinetic 
energy density functional Tg [p] . This gives the kinetic en- 
ergy of the ground state of a reference or model system of 
N noninteracting electrons in a self-consistent field cho- 
sen so that the ground state density equals p(r). The sec- 
ond term is the classical (Hartree) electron- electron po- 
tential energy (atomic units are used in this paper), and 
the last term is the exchange- correlation density func- 
tional Exc [p, ] which formally accounts for all remaining 
smaller corrections. 



If these functionals were known, then one could ob- 
tain the density p(r) from the variational principle (Euler 
equation) associated with eq 



p = VTAr; M) + T4//(r; [p]). 



(3) 



and then determine the total energy of the inhomoge- 
neous system from the energy functional E[p\. All other 
physical quantities related to the ground-state density 
could also be determined. Here 



(4) 



is the KS kinetic potential |9| , formally defined as a func- 
tional derivative of [p] , p is the chemical potential (the 
Lagrange multiplier associated with the normalization 
condition J p{r)dr — N), and Veff{r; [p]) is an effective 
one-body potential for the KS partitioning, defined by 



p(r') 



r — r' 



■dr' 



Sp{r) ■ 



(5) 
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Through the efforts of many workers, we now have 
rather accurate expressions for the exchange-correlation 
potential 14c (r; [p]) = SExc[p]/Sp{r) in eq and the 
main problem in solving the Euler equation is in find- 
ing an accurate approximation for the kinetic potential as 
a functional of the density 1, 10]. The most straightfor- 
ward way to proceed would be to develop approximations 
for Ts[p] itself, and obtain VT^{r] [p]) by functional differ- 
entiation, but this has proved to be a very challenging 
problem. As discussed below, simple local density ap- 
proximations for Ts[p] based on the Thomas- Fermi (TF) 
theory are not very accurate, and make several qualita- 
tively incorrect predictions. 

To bypass such problems, KS suggested a different way 
to compute the numerical value of Ts exactly, not directly 
from the density itself, but by introducing a set of N 
orbitals satisfying the N coupled KS equations that de- 
scribe the model system 0, |g . Using these orbitals one 
can determine both Ts and the total density p(r). 
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This accurate treatment of the kinetic energy played 
a central role in the development of DFT as a quanti- 
tative method [3. The KS method has yielded very ac- 
curate ground-state densities and energies for many sys- 
tems using currently available -Exc[p] • In many cases even 
a simple local density approximation for Exc[p\ proves 
adequate, and we will use this in the results reported 
here. However, the numerical cost of self-consistently de- 
termining N orbitals rapidly increases for large N . An 
accurate treatment of the kinetic energy as well as the 
potential energy contributions in terms of the electron 
density only, as originally envisioned in the work of HK, 
would certainly be desirable. 

To that end, there has been renewed interest in devel- 
oping "orbital-free" (OF) approximations for the kinetic 
energy density functional T5[p], with some notable suc- 
cesses in describing extended systems 0, 0, 0, 0, 0, 
[rH - The main advance has been to introduce nonlocal 
density functionals that reproduce known exact results 
for the linear response of the density of the model system 
to small perturbations. Similar ideas led to the develop- 
ment of nonlocal weighted density functional theories for 
classical nonuniform fluids, again with some notable suc- 
cesses However, in both cases it has proved difficult 
to gain a physical understanding of how best to incor- 
porate essential nonlocal effects directly into the kinetic 
energy (or classically, the free energy) functional through 
various weighting kernels and to understand what errors 
in the density will likely arise from given approximations 
to these functionals. 

We pursue here an intermediate strategy between the 
ultimately desirable goal of determining Ts[p] as a func- 
tional of the density and the KS method, which avoids 
direct use of the density and uses orbitals to calculate 
numerically both Ts and the density. We focus instead 
on approximating the kinetic potential Vt^{y\[p]) as a 
functional of the density. The connection to the density 
response is more direct and physically suggestive and we 
believe that nonlocal effects in VT^{'fc\ {p\) could be less 
problematic than in the full TsM- 



II. NONLOCAL EFFECTS IN T^p] AND Ft, (r; [p]) 

Ts[p\ can be formally related to Vt^{v; [p]) by a (func- 
tional) integration over the density changes in all space, 
as the density is changed from a uniform to the final 
nonuniform value 0, 0| . Because of this spatial inte- 
gration, Ts[p] is a more nonlocal functional of the density 
than is ^^-^(r; [p])- Other detailed theoretical arguments 
suggesting that Vt, is more local have recently appeared 
[l9| . A second consequence of this integration is that the 
numerical value of Ts [p] depends mainly on smoothed or 
spatially averaged properties of p(r). Thus the inverse 
problem of accurately determining p(r) from eq ^ at 
all values of r using a Vr, (r; [p]) derived from an approx- 
imate Ts[p\ in eq 10} can be poorly conditioned: very 
large errors in the density at particular values of r can 



arise from an approximate Ts[p] whose numerical value 
differs by only a relatively small amount from the correct 
value. For example, as discussed below, a self-consistent 
solution of eq ||2J) using a Vt^ (r; [p]) arising from the local 
Thomas- Fermi approximation to Tg [p] and a 14/ / (r) ap- 
propriate for an Ar atom predicts an infinite density at 
the nucleus, although the numerical value for Ts is only 
about 30% above the correct value. 

But there is another way to make progress. If we can 
find a reasonably accurate way to calculate p(r) directly 
from eq Q for a variety of T4a:t(r), then we can integrate 
the induced density changes as the field is "turned on" 
to obtain the total energy E and Tg. See eq (|16f) . Now 
the smoothing property of the integration helps make the 
results for E less sensitive to whatever errors remain in 
p(r) and we would expect good results [iTj . 

However, in general the integration must be carried out 
numerically over a range of density changes and simple 
analytical approximations to Vt^'^', [p\) do not automati- 
cally translate into equivalently simple analytical expres- 
sions for Ts [p] . If accurate direct approximations for [p] 
could be found the corresponding determination of E and 
Ts would require only the (final) equilibrium density p(r) 
and would certainly be simpler. It remains to be seen 
whether the benefits of focusing first on VTX^'^ \p\) '^iU 
outweigh this cost. In the interim our approach still pre- 
serves many of the potential advantages of OF methods 
over the more complicated KS procedure of determin- 
ing Ts numerically using orbitals, and we believe that 
insights gained from its use may lead to better direct 
approximations for Tg [p] . 

Similar ideas has been successfully applied to classical 
nonuniform fluids 0, 0, |23, 0| , where the density re- 
sponse and free energy change arising from a wide variety 
of external fields Vext (r) have been accurately determined 
from the classical analogue of eq (O . The terms in that 
equation can be given physically suggestive interpreta- 
tions and the results are comparable to those arising from 
the best versions of classical weighted density functional 
theories. In the simplest version of this theory, appropri- 
ate for very slowly varying effective fields, one arrives at 
a (local) differential equation involving only the density 
and a linear second derivative term, as in the classical 
van der Waals theory of the liquid- vapor interface [Tsji . 
Nonlocal effects arising from more rapidly varying fields 
can be taken into account directly in the density response 
by usin g a.n op timized version of linear response theory 

mrnmm 

In this paper we consider only local approximations to 
VtX^'i [p]) that at each r involve the local density p(r) 
and its lowest order gradients, motivated by analogous 
approximations made for Tg [p] in the Thomas-Fermi the- 
ory and the von Weizsacker modifications and by the 
equivalent theory for classical fluids. In future work we 
plan to examine nonlocal corrections suggested by an op- 
timized version of linear response theory and thus make 
contact with current work on OF methods and earlier 
related work by Pratt and coworkers p3 |. 
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III. THOMAS-FERMI AND VON WEIZSACKER 
APPROXIMATIONS FOR Ts[p] 

The first direct approximation for the kinetic energy 
functional was found in the Thomas-Fermi (TF) model 
mill, where 

Ttf[p\^Cp j pi{r)dv. (6) 

Here Cp = ^(37r^)i. The TF kinetic energy functional, 
Ttf[p], is derived by local use of a uniform free electron 
gas model and is known to be exact in the limit of an 
infinite number of electrons [2^ . 

However, for an atomic system, the number of elec- 
trons N is finite and the effective potential is very rapidly 
varying near the nucleus. As a result, there are many de- 
ficiencies in the TF model for such applications. The 
total energies predicted for atomic systems give at best 
an order of magnitude estimate, and cannot be relied on 
for quantitative calculations. More serious errors can be 
seen in the density response at particular values of r as 
determined from the Euler equation Q. As mentioned 
above, the TF model predicts an infinite electron density 
at the nucleus. In addition, the density does not decay 
exponentially in the classically forbidden region (the tail 
region), and it does not show any shell structure. Fur- 
thermore, the TF model predicts no atomic bonding to 
form molecules or solids 26, 27j. 

Though the TF model has many shortcomings, it is a 
natural local approximation and gives a simple analytical 
expression for the kinetic energy density functional and 
the kinetic potential. Moreover, it is known that greatly 
improved results can be obtained for atomic and molecu- 
lar ground state energies [2^ |2^ |^ |^ if better quality 
densities are used in Ttf[p\ instead of the self-consistent 
variational densities arising from a solution of eq For 
example, if the Hartree-Fock (HF) density [s^l for Ar is 
used in Ttf[p], the kinetic energy is only 7% below the 
correct result. This suggests that improvements to the 
TF model should focus directly on the Euler equation (O 
and the density response it predicts. 

In this paper we show that the Euler equation ||2Jl can 
give surprisingly accurate results for atomic densities and 
energies if the kinetic potential derived from the original 
TF model is augmented by a simple gradient term that 
helps satisfy important boundary conditions both near 
and far away from the nucleus. The specific form of this 
correction was inspired by the von Weizsacker modifica- 
tion of the TF kinetic energy density functional, which 
we now describe. 

By considering modified plane waves of the form (1 + 
a ■ r)e*'^ '", with a a constant vector and k the local wave 
vector, von Weizsacker |^ derived a correction 

to the TF kinetic energy functional. A generalized total 
kinetic energy functional incorporating a correction of 



this kind can be written as 

Ttfxw [p] = Ttf [p] + ATiv [p] , (8) 

where A is a parameter. A = 1 corresponds to the original 
von Weizsacker correction, and its self-consistent use in 
eq ^ for atoms not only yields a finite electron density 
at the nucleus with the exact cusp condition satisfied, 
but also gives the correct exponential decay of the density 
far from the nucleus [SSLISq for the exact p. We refer to 
this as the TFW model. Disappointingly, however, the 
kinetic energy given by this approach is not much better 
than that of the original TF model. For example, the 
results for Ar are 29% below the correct value. 

Other values of A have also been examined. All TFAW 
models predict a finite electron density at the nucleus, 
and exponential decay far from nucleus but only A = 1 
satisfies the exact cusp condition and gives the correct 
exponential decay constant. A = 1/9 has been shown 
to give the correct second order gradient expansion cor- 
rection to the TF functional and to reproduce the ex- 
act linear response of a homogeneous electron gas under 
a long wavelength perturbation but energy predic- 
tions for atoms again are inaccurate. Tomishima and 
Yonei found empirically that rather accurate atomic 
ground state energies for a wide range of Z can be ob- 
tained by using A = 1/5. However the inaccurate descrip- 
tion of the electron density near and far from nucleus and 
the absence of shell structure are still present, and part of 
this good agreement probably comes from a cancellation 
of errors. 

These models use simple gradient corrections sug- 
gested by the von Weizsacker theory to try to improve the 
TF kinetic energy functional directly. However nonlocal 
effects certainly must be important in Ts[p], as is recog- 
nized in recent OF methods, and these must be taken 
into account before this approach is likely to be quanti- 
tatively useful. 

IV. MODIFIED EULER EQUATIONS FOR THE 
DENSITY RESPONSE 

As argued above, it may be more profitable to exam- 
ine simple gradient corrections directly in the more local 
Euler equation. The Euler equation arising from use of 
the TFAW expression for Ts [p] is 

5^ 2^ ^ A |Vp(r)|2 A V2p(r) , , , . 

There is a local term in the density from the TF model 
and two terms involving density gradients. Both terms 
must be present with specified coefficients if we insist on 
using simple local expressions of the von Weizsacker form 
directly in [p] with no compensating terms from nonlo- 
cal corrections. However the first nonlinear gradient term 
does not have a simple physical interpretation and does 
not appear in the analogous classical van der Waals equa- 
tion for the liquid- vapor interface [3| ■ A linear gradient 
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correction like the second term would be expected on 
quite general grounds. Moreover we find that the second 
term alone is responsible for producing a finite density 
at an atomic nucleus with Vext{ic) = Vextif) = —Z/r, 
effectively building in the uncertainty principle, and of 
course this term alone contributes to the linear response 
of the density to weak fields. Finally this term yields 
exponential decay far from the nucleus. 

Thus we propose as the simplest correction to the TF 
kinetic potential a modified Euler equation where only 
the linear gradient term appears, with a coefficient chosen 
to best describe the density response. To avoid confusion 
in notation, A is changed to a in our new Euler equation: 



5 „ 2 a V2p(r) 
3 4 p(r) 



+ K//(r). (10) 



Note that we use a Ve//(r) that includes the effects of 
exchange and correlation in all the calculations reported 
herein. Thus Eq. (|10|l does not reduce to the classical 
TF model when a = 0. 

We first examine the predictions of this model for 
atoms. It is easy to show that the density at the nu- 
cleus is finite for any choice of a and satisfies 



p{0) 



2Z 



a 



(11) 



The exact cusp condition |3J| for the logarithmic deriva- 
tive of the density at the origin can be satisfied if a = 1, 
as in the TFW model. However this choice does not 
guarantee correct values for p(0) and p'(0) separately, 
and indeed the TFW model gives very poor results for 
these terms. For example, for Ar p{Q) is 72% below the 
HE value. Moreover the nonlinear gradient term in the 
TFW model does not seem to improve the density at in- 
termediate values of r when compared to the results of 
our model. 

Equation H1U|) also produces exponential decay of the 
density at large r, where 



p{r) w exp[— (— 4/i/Q;)2r] 



When a = 1/2, the correct exponential decay constant 
for the exact p, is predicted. 

Using this one parameter model we cannot satisfy both 
the exact cusp condition and asymptotic decay condition 
with the same a, and it is clear that no oscillatory shell 
structure will be generated by the simple gradient term. 
(As is the case for the classical van der Waals theory of 
interfaces, nonlocal corrections are required to get oscil- 
latory structure |l8j|.') But we can hope that a proper 
choice of a could give a smoothed electron density that 
reproduces important overall features and gives improved 
atomic energies when compared to the TF model. Since 
any a produces a finite density at the nucleus and thus 
corrects a major shortcoming of the TF model, we choose 
a = 1/2 to give the correct exponential decay of the den- 
sity far from the nucleus. This is the region that should 



be most important in chemical bonding. In the follow- 
ing we refer to eq ffTU\i with a = 1/2 as the modified 
Thomas-Fermi (MTF) model. We show below that the 
MTF model gives atomic energies comparable to those of 
the best TFAW model with A = 1/5 along with a much 
better description of the density distribution. 



V. RESULTS FOR THE MTF MODEL 

A. Linear response function 

We first examine the predictions of the various models 
for the linear response of the density of the model system 
to a weak perturbing potential. In Fourier space this can 
be written as 



6p{k)^x{^)SV{k). 



(13) 



Here x(k) is the Fourier transform of the linear response 
function for a uniform system of noninteracting electrons, 
which can be derived exactly from first-order perturba- 
tion theory ^| : 



/, N f^f A 1 — 9^ 1 I 1 + 9 l^ 

X k = Tj + ^^ln - 

TT^ 2 Aq 1 — 9 



(14) 



where kf is the Fermi wave vector and q = k/2kf is a 
dimensionless momentum. 

The corresponding response function for the TFAW 
model, eq ©, is given by ^3 



X(k) 



TFXW 



7r2'l + 3A(?' 



r), 



(15) 



as can be seen from eq (O for an infinitesimal Ve//. It 
is easy to see that the nonlinear gradient term in Q 
does not contribute to the linear response function. As a 
result, the linear response function for the MTF model, 
derived from eq (|10|l with a = 1/2, is the same as the 
TFiW model. 

As shown in Figure 1, the TF^W model is exact for 
long- wavelength perturbations, but becomes worse in the 
high-fc region, while the TFW model is exact for short- 
wavelength perturbations, but fails in the \ow-k regions. 
The MTF model gives a qualitatively reasonable aver- 
age description of the exact linear response function, es- 
pecially in the important region near the singularity at 
q = 1. Wang et al. J^] have suggested that reproducing 
the singularity and the overall form of the linear response 
function are important for obtaining the correct physics 
from OF methods, particularly for producing the shell 
structure of atoms. While no simple gradient type model 
can describe the singularity exactly (or give shell struc- 
ture), the MTF model captures the average behavior of 
the linear response function and shows better accuracy 
in the physically significant region near q = 1. 
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FIG. 1: Linear response functions of a uniform system of 
independent fermions as given by the various models. 



B. Atoms 



TABLE I: Electron density at the nucleus p(0) using the 
various TF-type models, the KS method, the HF method and 
the MTF model. 





MTF 




TF^W 


TFW 


KS 


H 


0.5820 


0.3183 


2.390 


0.1029 


0.2728 


He 


5.208 


3.596 


19.29 


0.9515 


3.525 


Ne 


788.2 


619.9 


2595 


169.6 


614.5 


Ar 


4811 


3840 


15482 


1093 


3819 


Kr 


40064 


32236 


126491 


9579 


32146 


Xe 


137631 


112219 


431005 


33710 


111005 



"the HF data is taken from Clementi et al l32| . 



TABLE IL (r-i) using the TFfW model, the KS method, 
and the MTF model. 





MTF TF^W 


KS 


H 


1.000 


1.249 


0.9208 


He 


2.907 


3.261 


3.312 


Ne 


30.56 


30.75 


31.00 


Ar 


69.93 


69.63 


69.62 


Kr 


183.5 


181.7 


182.7 



Xe 321.2 317.8 317.8 



We carried out benchmark calculation on the hydro- 
gen atom and the rare gas atoms based on our MTF 
model. Our results are compared with the TP model, 
and its von Weizsacker-type gradient corrections. We 
examined A = ^ , ^, and 1 but report results only for 
A = i and 1, and for the KS and HF methods. The lo- 
cal density approximation (LDA) jU Elil El for the 
exchange-correlation functional is used for all the models. 
The standard method for solving Euler equations derived 
from the TP AW model 45] was implemented to solve the 
Euler equation for the MTF model and the KS theory. 
Our TPiW and KS results are in very good agreement 
with previous calculations 46] using a slightly different 
LDA parametrization The code uses the finite dif- 
ference method with the Gauss-Chebyshev (of the second 
kind) radial quadrature, proposed by Becke et at, con- 
sisting of 1000 points, and Becke's algorithms for solving 
Poisson's equation for the Hartree potential j48|]. 

As shown in Table the electron density at the nu- 
cleus, p{0), calculated by the MTF model is close to the 
HF results. As discussed above, the TP model predicts 
an infinite value for the electron density, while the TF-iW 
models overestimates the density by a factor of 4. The 
TFW model considerably underestimates the density, de- 
spite satisfying the exact cusp condition. It is known [49l ] 
that p(0) is related to the density moments (r~^) and 
(r~^). Thus the MTF model would be expected to give 
good results for these expectation values, as is shown in 
Table |n] and IIIII respectively. 

In Table we compare the chemical potential fi for 
the MTF theory, calculated by a self-consistent solution 
of eq H10(l . with that given by the KS equation and find 
very good agreement. This means that the asymptotic 
behavior of the density predicted by the MTF model is 
close to the prediction of the KS theory, since the MTF 



model would give exact results for the exact fi. Since the 
accuracy of the density in the tail is an important indica- 
tor of chemical bonding l] , we believe that this desired 
behavior could lead to chemical bonding in molecules. 

The Euler equation for the MTF model is not de- 
rived by a functional derivative of an energy functional. 
As a result we cannot immediately obtain the energy 
once we have determined the density, as would be pos- 
sible in OF methods that give direct approximations for 
TJp]- Instead, we can use a coupling parameter approach 
'22', Is^, 'FT, 52, "53] to calculate the energy change as the 
external potential is turned on: 

E = Ex=o+fJ. J [pA=i(r)-PA=o(r)]«ir+ J '^^ J P\{r)VextdT 

(16) 

Here the external potential Vext is linearly scaled by a 
coupling parameter A (not to be confused with the A 
parameter in the TFAVF methods) and the corresponding 
density is pA- E is the total energy of the system, and 
E\^o is the total energy at the same n in the absence of 



TABLE IIL (r^^) using the TFiW model, the KS method, 
the HF method, and the MTF model. 





MTF TF^W 


KS 


HF" 


H 


2.394 


4.720 


1.749 




He 


11.50 


20.06 


11.71 


11.99 


Ne 


394.7 


599.7 


411.9 


414.9 


Ar 


1390 


2056 


1459 


1465 


Kr 


6025 


8705 


6318 




Xe 


14101 


20148 


14799 


14818 



"the HF data is taken from Porras et al l57l . 
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TABLE IV: Atomic chemical potential fi using the various 
TF-type models, the KS method, and the MTF model. 





MTF 


TF 


TF^;W 


TFW 


KS 


H 


-0.2825 


-0.08039 


-0.09549 


-0.09773 


-0.2337 


He 


-0.3822 


-0.06290 


-0.1013 


-0.1393 


-0.5702 


Ne 


-0.4400 


-0.04959 


-0.1093 


-0.2188 


-0.4978 


Ar 


-0.4304 


-0.05034 


-0.1112 


-0.2346 


-0.3823 


Kr 


-0.4224 


-0.05017 


-0.1130 


-0.2485 


-0.3464 


Xe 


-0.4192 


-0.05034 


-0.1138 


-0.2552 


-0.3100 



TABLE V; Atomic energy E using the various TF-type mod- 
els, the KS method, and the MTF model. 





MTF 


TF 


TF^W 


TFW 


KS 


H 


-0.6092 


-1.076 


-0.6085 


-0.2927 


-0.4459 


He 


-2.902 


-4.756 


-2.917 


-1.559 


-2.834 


Ne 


-129.5 


-176.6 


-129.5 


-86.40 


-128.2 


Ar 


-529.0 


-681.1 


-526.2 


-375.5 


-525.9 


Kr 


-2774 


-3378 


-2747 


-2099 


-2750 


Xe 


-7293 


-8643 


-7214 


-5701 


-7229 



the external potential. 

If exact values for p\ were used, eq (|16|) would give ex- 
act results and a variety of different integration pathways 
(e.g., nonlinear scaling of Vext) could be used. Since we 
make approximations there will be errors in our predicted 
E values and different pathways could (incorrectly) give 
different results. However, our experience for the analo- 
gous classical problem suggests that these errors are small 
if the density is reasonably accurate ^3 1 as seems to be 
the case here. 

We note that p\=o and thus i?A=o vanish for the uni- 
form electron gas because of the fixed negative value of 
the chemical potential. Therefore, E, the total energy, 
can be straightforwardly computed from eq H16|) by nu- 
merical integration over a series of A values. For the 
benchmark calculations reported here for the total en- 
ergy, 1000 points are used between A = and A = 1, 
though reasonable results can be obtained with about 10 
points. The kinetic energy can then be found by subtract- 
ing the potential energy (calculated from the potential 
energy density functional) from the total energy. This 
way of obtaining kinetic energy is equivalent to integrat- 
ing the kinetic potential directly as discussed above, but 
we believe it is numerically more accurate. As shown in 
Table El the total energy of the MTF model is compara- 
ble to that of the TFiW model, the best TFAW model 
for the total atomic energy. 

In Figure 2-5, we compare the radial density distribu- 
tion r"^ p{r) of the MTF model to that predicted by other 
theories. Although the MTF model shows no shell struc- 
ture when compared to the essentially exact KS theory, it 
gives a nice averaging of the electron density, with slightly 
more emphasis on the principle peak than in the TFAW 
models, along with significantly better results both near 
and far from the nucleus. The fact that a smooth density 
can give such accurate energies illustrates the point that 
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FIG. 2: Radial density r^p of the Neon atom with the various 
TF-type models, the KS method, and the MTF model. 
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FIG. 3: Same as in Figure 2 but for the Argon atom 

the integration in eq (|16|) renders it less sensitive to small 
errors in the density. 

C. Diatomic molecules 

Gordon and Kim (GK) suggested a simple model for 
calculating the intermolecular potential of closed-shell 
molecules [s^. By assuming the molecular density of 
a diatomic molecule can be expressed as the sum of the 
two separate atomic densities, reasonable bonding poten- 
tials for closed-shell molecules can be obtained. Instead 
of using the inaccurate TF atomic densities, GK used the 
atomic HF densities and inserted these into the original 
TF kinetic energy functional, together with the classical 
electron-electron and nuclear-electron potentials, and a 
local density approximation for the exchange-correlation 
functionals. 

We followed GK's simple approach, using atomic densi- 
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FIG. 4: Same as in Figure 2 but for the Krypton atom FIG. 6; Ar-Ar interaction potential via the Gordon-Kim ap- 

proach. 




FIG. 5: Same as in Figure 2 but for the Xenon atom 



ties for the Ar atom given by the various models discussed 
here. Becke's atomic partition, the Gauss-Chebyshev (of 
the second kind) radial quadrature, and the Clenshaw- 
Curtis angular quadrature with 1000 radial points and 
500 ang ular points are used in the numerical calculations 

We found that the KS bonding potential is very simi- 
lar to the HF result, as is expected due to the high qual- 
ity of its atomic density. As we can see in Figure 6, the 
strength of the bonding of the TFAW models increases as 
A increases. However, the bonding energy is consistently 
overestimated, even for the smallest value of A. The MTF 
model predicts a reasonable bonding curve, with both its 
strength and bond length close to the HF and KS meth- 
ods. Note that even the TF density can give a quali- 
tatively accurate bonding curve using the GK method, 
though a self-consistent solution of the molecular density 
using the TF model is known to give no bonding [26] . 



VI. POSSIBLE GENERALIZATIONS 

We used heuristic arguments to suggest the form of the 
MTF equation (^UJ- One might hope that a more sys- 
tematic search using generalized forms of this kind could 
yield even better results. Since the MTF model has only 
one parameter, a, we considered one such generalization 
of the gradient term in H10|l that allows for the intro- 
duction of a second parameter, (3, and that interpolates 
between the MTF and the TFAW models: 



4 



f3iP - 1) 



|Vp(r) 



V2p(r) 



p(r)2 ' p{r) 

(17) 

For a = 1/2 and /3 = 1, this reduces to the MTF 
form, while the TFAW models correspond to {a, (3} = 
{2A, 1/2}, as can be seen from eq 

Using this generalized form in eq H1U|I for atomic sys- 
tems, the density at the nucleus is finite and satisfies a 
modified cusp condition 



P(0) 



2Z 



(18) 



When af3 — 1, the exact cusp condition is satisfied. This 
form also produces an exponentially decaying density in 
the tail, satisfying 



p{r) w exp 



2a/?2 



(19) 



When 2a/3^ — 1 the correct asymptotic decay constant 
is found for the exact /i. 

However only when {a, /?} — {2, 1/2} (i.e., the original 
von Weizsacker model with A = 1) can both conditions 
be satisfied simultaneously. Since all values give a finite 
density at the origin, we also considered the one parame- 
ter family of models with 2af3'^ = 1 that give the correct 
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exponential decay constant. However, the best overall 
results for the atomic energy and density were still found 
using the MTF model with {a, (3} = {1/2, 1}. Evidently 
the nonlinear gradient term in eq (|17|1 that is generated 
when (3^1 does not improve the simple MTF model 
and other forms should be examined. 



VII. FINAL REMARKS 

In summary, we have argued that it may be profitable 
to focus first on the Euler equation determining the elec- 
tron density when correcting the TF model and more 
generally in developing OF approximations in DFT. Non- 
local effects may be less severe, as is suggested by the 
generally good results found from use of the simple gra- 
dient correction discussed here. In future work we hope 
to take account of nonlocal corrections to the density re- 
sponse by using an optimized version of linear response 
theory 

As presently implemented in eq (|16|l . a coupling pa- 
rameter integration is required to determine the total 
energy directly from the density response. This is nu- 
merically more costly than would be the case if accurate, 
variationally stable, direct OF approximations for Ts[p\ 
can be found that yield good results for p(r) for all r when 
the associated VT^{r; [p]) is used in the Euler equation. 

A less ambitious goal that may prove almost as use- 
ful is to develop approximate energy functionals that. 



while not variationally accurate enough to use directly 
in the Euler equation, still give good numerical results 
for the kinetic and total energies of particular classes of 
systems if reasonably accurate densities are supplied as 
input. The use of such functionals in theories like ours 
would eliminate the need for the coupling parameter in- 
tegration and simplify the calculation of the total energy. 
As the simplest such example, if the MTF density is di- 
rectly inserted into the TF kinetic energy functional, the 
deviation of the kinetic energy from the MTF value is 
less than 10% (except for the lightest atoms H (20%) 
and He (17%)). We found empirically that the TF^W 
functional gives results to within about 1% for the heav- 
ier atoms, though there is no real justification for such 
a choice and a more systematic procedure is called for. 
We are currently investigating whether these ideas can 
be further developed and usefully applied to molecules 
and extended systems. 
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